#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <Source_CPgrid/mesh.hpp>

#include "PickingSpots.h"
#include "Source_Dense/DenseSolver.h"
#include "RBF_function.h"
#include "RBF_interpolation.h"

namespace rbf_interp {

    //获取最大和最小的属性值
    vector<double> getExtremum(const vector<double>& values) {
        double min_value = *min_element(values.begin(), values.end());
        double max_value = *max_element(values.begin(), values.end());
        return { min_value, max_value };
    }
    // 归一化函数
    vector<double> normalize(const vector<double>& values) {
        double min_value = getExtremum(values)[0];
        double max_value = getExtremum(values)[1];
        vector<double> normalized_values;
        for (double value : values) {
            normalized_values.push_back((value - min_value) / (max_value - min_value));
        }
        return normalized_values;

    }
    // 归一化逆操作函数
    vector<double> denormalize(const vector<double>& normalized_values, double min_value, double max_value) {
        vector<double> denormalized_values;

        for (double value : normalized_values) {
            denormalized_values.push_back(value * (max_value - min_value) + min_value);
        }
        return denormalized_values;

    }
    double denormalize_number(double normalized_values, double min_value, double max_value) {
        double denormalized_values = normalized_values * (max_value - min_value) + min_value;
        return denormalized_values;
    }

    //计算向量叉积
    vector<double> crossProduct(const vector<double>& v1, const vector<double>& v2) {
        return {
            v1[1] * v2[2] - v1[2] * v2[1],
            v1[2] * v2[0] - v1[0] * v2[2],
            v1[0] * v2[1] - v1[1] * v2[0]
        };
    }

    // 计算向量点积
    double dotProduct(const vector<double>& v1, const vector<double>& v2) {
        return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
    }

    // 计算四面体体积
    double tetrahedronVolume(const vector<double>& a, const vector<double>& b, const vector<double>& c, const vector<double>& d) {
        vector<double> ab = { b[0] - a[0], b[1] - a[1], b[2] - a[2] };
        vector<double> ac = { c[0] - a[0], c[1] - a[1], c[2] - a[2] };
        vector<double> ad = { d[0] - a[0], d[1] - a[1], d[2] - a[2] };
        vector<double> cross = crossProduct(ab, ac);
        double volume = abs(dotProduct(ad, cross)) / 6.0;
        return volume;
    }

    // 计算六面体总体积
    double hexahedronVolume(const vector<vector<double>>& points) {
        // 计算六面体总体积
        //底面顶点坐标: points[0], points[1], points[2], points[3]
        //顶面顶点坐标: points[4], points[5], points[6], points[7]
        double V1 = tetrahedronVolume(points[0], points[1], points[2], points[4]);
        double V2 = tetrahedronVolume(points[1], points[2], points[3], points[4]);
        double V3 = tetrahedronVolume(points[2], points[3], points[0], points[4]);
        double V4 = tetrahedronVolume(points[3], points[0], points[1], points[4]);
        double V5 = tetrahedronVolume(points[4], points[5], points[6], points[0]);
        double V6 = tetrahedronVolume(points[5], points[6], points[7], points[0]);
        double V = V1 + V2 + V3 + V4 + V5 + V6;
        return V;
    }


    //计算守恒系数的值
    vector<double> calculatedConservationCoefficient(
        double total_volume,
        double total_rockden,
        double total_poro,
        const vector<double>& volume,
        const vector<double>& rockden,
        const vector<double>& poro) {
        //总质量
        double m_total = total_rockden * total_volume;
        //总的孔隙体积
        double v_total = total_poro * total_volume;
        //求和质量
        double m_count = 0;
        //求和体积
        double v_count = 0;
        //
        double v_sum = 0;
        double flag = 0;
        for (size_t i = 0; i < volume.size(); ++i) {
            m_count = m_count + (rockden[i] * volume[i]);
            v_count = v_count + (poro[i]* volume[i]);
            v_sum += volume[i];
            flag++;
        }
        double coefficient_m = m_total / m_count;
        double coefficient_v = v_total / v_count;
        return {coefficient_m, coefficient_v};
    }


    double rbfGaussian(double distance) {
        // RBF形状参数
        double epsilon = 0.01;
        return exp(-(epsilon * distance) * (epsilon * distance));
    }

    // 半方差函数，这里使用的是线性模型
    double semivariance(double h, double nugget, double sill, double range) {
        if (h <= range) {
            return nugget + sill * (1.5 * (h / range) - 0.5 * pow(h / range, 3));
        } else {
            return nugget + sill;
        }
    }

    double euclideanDistance(
        const vector<double>& firstPointCoord,
        const vector<double>& secondPointCoord) {

        double sum = 0.0;
        for (size_t i = 0; i < firstPointCoord.size(); ++i) {
            sum += (firstPointCoord[i] - secondPointCoord[i]) * (firstPointCoord[i] - secondPointCoord[i]);
        }
        return sqrt(sum);
    }


    void denormalize(vector<double>& normalized_values, const vector<double>& values)
    {
        //获取values值向量最大最小值
        auto values_range = getExtremum(values);
        double min_value = values_range[0];
        double max_value = values_range[1];

        //对normalized_values进行逆归一化
        for (int i = 0; i < normalized_values.size(); ++i)
        {
            normalized_values[i] = normalized_values[i] * (max_value - min_value) + min_value;
        }

    }


    bool invertMatrix(
        vector<vector<double>>& matrix,
        size_t n,
        vector<vector<double>>& inverseMatrix) {

        vector<vector<double>> tempMatrix = matrix;
        inverseMatrix = vector<vector<double>>(n, vector<double>(n, 0.0));

        // 初始化逆矩阵为单位矩阵
        for (size_t i = 0; i < n; ++i) {
            inverseMatrix[i][i] = 1.0;
        }

        // 高斯-约当消元法
        for (size_t i = 0; i < n; ++i) {
            // 寻找主元
            double maxElement = abs(tempMatrix[i][i]);
            size_t maxRow = i;
            for (size_t k = i + 1; k < n; ++k) {
                if (abs(tempMatrix[k][i]) > maxElement) {
                    maxElement = abs(tempMatrix[k][i]);
                    maxRow = k;
                }
            }

            // 如果主元为零，则矩阵不可逆
            if (maxElement == 0) {
                return false;
            }

            // 交换行
            swap(tempMatrix[i], tempMatrix[maxRow]);
            swap(inverseMatrix[i], inverseMatrix[maxRow]);

            // 归一化当前行
            for (size_t j = i + 1; j < n; ++j) {
                tempMatrix[i][j] /= tempMatrix[i][i];
            }
            for (size_t j = 0; j < n; ++j) {
                inverseMatrix[i][j] /= tempMatrix[i][i];
            }

            // 消元
            for (size_t k = 0; k < n; ++k) {
                if (k != i) {
                    double factor = tempMatrix[k][i];
                    for (size_t j = i; j < n; ++j) {
                        tempMatrix[k][j] -= factor * tempMatrix[i][j];
                    }
                    for (size_t j = 0; j < n; ++j) {
                        inverseMatrix[k][j] -= factor * inverseMatrix[i][j];
                    }
                }
            }
        }

        return true;
    }

    vector<double> vectorMatrixMultiply(
        const vector<double>& valueVector,
        const vector<vector<double>>& inverseMatrix) {

        //得到逆矩阵的行rows
        size_t rows = inverseMatrix.size();
        //得到逆矩阵的列
        size_t cols = inverseMatrix[0].size();
        //声明结果容器
        vector<double> result(rows, 0.0);
        //计算矩阵和向量相乘
        for (size_t i = 0; i < rows; ++i) {
            for (size_t j = 0; j < cols; ++j) {
                result[i] += valueVector[j] * inverseMatrix[j][i];
            }
        }
        return result;
    }

    double rbfInterpolateValue(
        const vector<vector<double>>& roofCoood,
        const vector<double>& sValues,
        const vector<double>& targetPoint) {

        //size_t表示贡献点的维度
        size_t numPoints = roofCoood.size();
        //声明插值矩阵A
        vector<vector<double>> A(numPoints, vector<double>(numPoints, 0.0));
        //声明值向量b
        vector<double> b(numPoints, 0.0);

        //声明权重向量coefficients
        vector<double> coefficients(numPoints, 0.0);

        //声明插值矩阵A的逆矩阵inverseA
        vector<vector<double>> inverseA(numPoints, vector<double>(numPoints, 0.0));

        // 构建插值矩阵和右侧向量
        for (size_t i = 0; i < numPoints; ++i) {
            for (size_t j = 0; j < numPoints; ++j) {
                double r = euclideanDistance(roofCoood[i], roofCoood[j]);
                A[i][j] = rbfGaussian(r);
            }
            b[i] = sValues[i];

        }

        // 求解矩阵A的逆
        if (!invertMatrix(A, numPoints, inverseA)) {
            cerr << "Matrix A is singular and cannot be inverted." << endl;
            return {}; // 返回空向量表示失败
        }

        // 计算插值系数
        coefficients = vectorMatrixMultiply(b, inverseA);

        // 使用插值函数计算插值点的值
        double interpolatedValue = 0.0;

        for (size_t i = 0; i < numPoints; ++i) {
            double r = euclideanDistance(targetPoint, roofCoood[i]);
            interpolatedValue += coefficients[i] * rbfGaussian(r);

        }
        //返回拟合的值（应力分量）
        return interpolatedValue;
    }

    vector<double> rbfInterpolateTwovalues(
        const vector<vector<double>>& roofCoood,
        const vector<vector<double>>& values,
        const vector<double>& targetPoint) {

        //size_t表示贡献点的维度
        size_t numPoints = roofCoood.size();
        //声明插值矩阵A
        vector<vector<double>> A(numPoints, vector<double>(numPoints, 0.0));
        //声明第一个值向量b1
        vector<double> b1(numPoints, 0.0);
        //声明第二个值向量b2
        vector<double> b2(numPoints, 0.0);
        //声明权重向量coefficients1、coefficients2
        vector<double> coefficients1(numPoints, 0.0);
        vector<double> coefficients2(numPoints, 0.0);
        //声明插值矩阵A的逆矩阵inverseA
        vector<vector<double>> inverseA(numPoints, vector<double>(numPoints, 0.0));

        // 构建插值矩阵和右侧向量
        for (size_t i = 0; i < numPoints; ++i) {
            for (size_t j = 0; j < numPoints; ++j) {
                double r = euclideanDistance(roofCoood[i], roofCoood[j]);
                A[i][j] = rbfGaussian(r);
            }
            b1[i] = values[i][0];
            b2[i] = values[i][1];
        }

        // 求解矩阵A的逆
        if (!invertMatrix(A, numPoints, inverseA)) {
            cerr << "Matrix A is singular and cannot be inverted." << endl;
            return {}; // 返回空向量表示失败
        }

        // 计算插值系数
        coefficients1 = vectorMatrixMultiply(b1, inverseA);
        coefficients2 = vectorMatrixMultiply(b2, inverseA);

        // 使用插值函数计算插值点的值
        double interpolatedValue1 = 0.0;
        double interpolatedValue2 = 0.0;
        for (size_t i = 0; i < numPoints; ++i) {
            double r = euclideanDistance(targetPoint, roofCoood[i]);
            interpolatedValue1 += coefficients1[i] * rbfGaussian(r);
            interpolatedValue2 += coefficients2[i] * rbfGaussian(r);
        }
        //返回拟合的值（杨氏模量和泊松比）
        return { interpolatedValue1, interpolatedValue2 };
    }

    /*
       RBF插值函数
       points：贡献点坐标
       values：贡献点对应的杨氏模量、泊松比、岩石密度、孔隙度
       queryPoint：目标点坐标
    */
    vector<double> rbfInterpolateFourValues(
        const vector<vector<double>>& points,
        const vector<vector<double>>& values,
        const vector<double>& queryPoint) {

        //size_t表示贡献点的维度
        size_t numPoints = points.size();
        //声明插值矩阵A
        vector<vector<double>> A(numPoints, vector<double>(numPoints, 0.0));
        //声明第一个值向量b1,b2,b3,b4
        vector<double> b1(numPoints, 0.0), b2(numPoints, 0.0), b3(numPoints, 0.0), b4(numPoints, 0.0);
        //声明权重向量coefficients1、coefficients2、coefficients3、cofficents4
        vector<double> coefficients1(numPoints, 0.0), coefficients2(numPoints, 0.0), coefficients3(numPoints, 0.0), coefficients4(numPoints, 0.0);
        //声明插值矩阵A的逆矩阵inverseA
        vector<vector<double>> inverseA(numPoints, vector<double>(numPoints, 0.0));

        // 构建插值矩阵和右侧向量
        for (size_t i = 0; i < numPoints; ++i) {
            for (size_t j = 0; j < numPoints; ++j) {
                double r = euclideanDistance(points[i], points[j]);
                A[i][j] = rbfGaussian(r);
            }
            //A[i][i] += 1.0;
            b1[i] = values[i][0];
            b2[i] = values[i][1];
            b3[i] = values[i][2];
            b4[i] = values[i][3];
        }
        // 求解矩阵A的逆
        if (!invertMatrix(A, numPoints, inverseA)) {
            cerr << "Matrix A is singular and cannot be inverted." << endl;
            return {}; // 返回空向量表示失败
        }

        // 计算插值系数
        coefficients1 = vectorMatrixMultiply(b1, inverseA);
        coefficients2 = vectorMatrixMultiply(b2, inverseA);
        coefficients3 = vectorMatrixMultiply(b3, inverseA);
        coefficients4 = vectorMatrixMultiply(b4, inverseA);
        // 使用插值函数计算插值点的值
        double interpolatedValue1 = 0.0;
        double interpolatedValue2 = 0.0;
        double interpolatedValue3 = 0.0;
        double interpolatedValue4 = 0.0;
        for (size_t i = 0; i < numPoints; ++i) {
            double r = euclideanDistance(queryPoint, points[i]);
            //r = r / maxDistance;
            interpolatedValue1 += coefficients1[i] * rbfGaussian(r);
            interpolatedValue2 += coefficients2[i] * rbfGaussian(r);
            interpolatedValue3 += coefficients3[i] * rbfGaussian(r);
            interpolatedValue4 += coefficients4[i] * rbfGaussian(r);
        }
        //返回拟合的值（杨氏模量、泊松比、岩石密度、孔隙度）
        return { interpolatedValue1, interpolatedValue2,interpolatedValue3,interpolatedValue4 };
    }


    vector<vector<double>> rbfExtrapolate(const vector<vector<double>>& roofCoord,
                                          const vector<vector<double>>& values,
                                          const vector<vector<double>>& targetPoints_ext, int thread)
    {
        // cout << "Extrapolation of extended regional attributes start." << endl;
        // cout << "The size of matrix A is " << roofCoord.size() << endl;
        // cout << "The size of targetPoints_ext is " << targetPoints_ext.size() << endl;
        //插值矩阵大小：matrix_size * matrix_size
        size_t matrix_size = roofCoord.size();
        size_t value_types = values.size();
        size_t target_size = targetPoints_ext.size();
        //正则化系数
        double lambda_reg = 0.1;
        //插值矩阵A
        DenseMatrix A(matrix_size,matrix_size,0.0);
        //值向量B
        vector<DenseMatrix> B(value_types, DenseMatrix(matrix_size, 1, 0.0));
        vector<vector<double>> result(value_types,vector<double>(target_size,0.0));

        //定义RBF核函数
        RBF_function func_("Gaussian");
        // 构建插值矩阵A
        clock_t time_stt = clock();
#pragma omp parallel for schedule(static) num_threads(thread)
        for (int i = 0; i < matrix_size; ++i) {
            for (int j = 0; j < matrix_size; ++j) {
                double r = euclideanDistance(roofCoord[i], roofCoord[j]);
                // A.setValue(i, j,rbfGaussian(r));
                A.setValue(i, j, func_.weights(r));
            }
            // 添加正则项
            A.setValue(i, i, A.getValue(i,i) + lambda_reg);

            // 构建向量b
            for (size_t j = 0; j < value_types; j++)
                B[j].setValue(i, 0, values[j][i]);
        }

        // cout << "time use in Construct the interpolation matrix A is: " << (clock() - time_stt) / (double)CLOCKS_PER_SEC << "s" << endl;

        // 求解线性方程组 A * X = B
        DenseSolver solver(A,thread);
        time_stt = clock();
        auto coefficients = solver.solve(B);
        // cout << "time use in Solve AX = B is: " << (clock() - time_stt) / (double)CLOCKS_PER_SEC << "s" << endl;

        //获取插值系数
        vector<vector<double>> coefficients_(value_types,vector<double>(matrix_size,0.0));
        for (int i = 0; i < value_types; i++)
        {
            coefficients_[i] = coefficients[i].GetData_c(0);
        }

        //计算RBF插值结果
        time_stt = clock();
#pragma omp parallel for schedule(static) num_threads(thread)
        for (int i = 0; i < target_size; ++i)
        {
            for (int j = 0; j < matrix_size; ++j)
            {
                auto r = euclideanDistance(targetPoints_ext[i], roofCoord[j]);
                auto weight = func_.weights(r);
                for (int k = 0; k < value_types; k++)
                {
                    result[k][i] += coefficients_[k][j] * weight;
                }

            }

        }
        // cout << "time use in get the interpolation results is       " << (clock() - time_stt) / (double)CLOCKS_PER_SEC << "s" << endl;
        return result;
    }

    void rbfInterpolate(CPgrid& ori_grid, finiteGrid& ext_grid, int thread_num)
    {

        vector<double> lambdaValues = { 0,1e-13,1e-12, 1e-11,1e-10,1e-9, 1e-8,1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1, 10,100};

        //插值矩阵构建选点
        PickingSpots pickingSpots(0.01);
        auto selection = pickingSpots.Picking(ori_grid.hexsize);
        //获取Reservior区域网格单元信息
        auto data = ori_grid.hexdata;
        //获取selection坐标集
        vector<vector<double>> roofCoord(selection.size(),vector<double>(3,0.0));
        for (int i = 0; i < selection.size(); i++)
        {
            roofCoord[i][0] = data[selection[i]].center.x;
            roofCoord[i][1] = data[selection[i]].center.y;
            roofCoord[i][2] = data[selection[i]].center.z;
        }
        //获取selection属性集
        vector<vector<double>> values(4,vector<double>(selection.size(),0.0));
        for (int i = 0; i < selection.size(); i++)
        {
            values[0][i] = data[selection[i]].poissonr;
            values[1][i] = data[selection[i]].youngmode;
            values[2][i] = data[selection[i]].rockden;
            values[3][i] = data[selection[i]].poro;
        }
        //获取储层区域网格单元数量
        auto cells_res = ext_grid.topo_i * ext_grid.topo_j * ext_grid.topo_k;
        //获取延拓区域网格单元数量
        auto cells_ext = ext_grid.i_withext * ext_grid.j_withext * ext_grid.k_withext - cells_res;
        auto& ext_data = ext_grid.hexdata;
        //计算剖分后的网格单元中心坐标
#pragma omp parallel for schedule(static) num_threads(thread_num)
        for (int i = 0; i < ext_grid.hexsize; ++i)
        {
            ext_data[i].cal_center(&ext_grid.pdata[0]);
        }
        //遍历延拓区域，获取延拓插值区域信息
        vector<vector<double>> targetPoints_res(cells_res, vector<double>(3));
        vector<int> target_index_res(cells_res);
        int cell_index = 0;
        int sum = 0;

        for (int k = 0; k < ext_grid.k_withext; ++k)
        {
            for (int i = 0; i < ext_grid.i_withext; ++i)
            {
                for (int j = 0; j < ext_grid.j_withext; ++j)
                {
                    if (i >= ext_grid.reservior_relative_i_start && i < ext_grid.reservior_relative_i_end &&
                        j >= ext_grid.reservior_relative_j_start && j < ext_grid.reservior_relative_j_end &&
                        k >= ext_grid.reservior_relative_k_start && k < ext_grid.reservior_relative_k_end)
                    {
                        cell_index = ext_grid.get_hex_id(i,j,k, ext_grid.i_withext, ext_grid.j_withext, ext_grid.k_withext);
                        targetPoints_res[sum][0] = ext_data[cell_index].center.x;
                        targetPoints_res[sum][1] = ext_data[cell_index].center.y;
                        targetPoints_res[sum][2] = ext_data[cell_index].center.z;
                        target_index_res[sum] = cell_index;
                        ++sum;
                    }
                }
            }
        }
        double max_rbf_poro = -10;

        //储层插值并守恒
#pragma omp parallel for schedule(static) num_threads(thread_num)
        for (int i = 0; i < ori_grid.hexsize; ++i) {
            //声明坐标集——用于存储原始网格插值贡献点
            vector<vector<double>> root_coord;
            vector<double> temp_coord;
            //声明原始网格坐标对应的值
            vector<vector<double>> root_values(4);
            //vector<vector<double>> root_values;
            vector<double> root_temp_values;

            //声明原始网格的ijk
            int root_i, root_j, root_k;
            //获得当前原始网格的ijk
            ori_grid.get_hex_ijk(i, root_i, root_j, root_k);

            //声名root_id,用于存储当前节点周围的id
            int root_id;
            vector<vector<double>> a_values;


            //获取周围网格单元的属性值
            for (int a = root_i - 1; a <= root_i + 1; ++a) {
                for (int b = root_j - 1; b <= root_j + 1; ++b) {
                    for (int c = root_k - 1; c <= root_k + 1; ++c) {
                        if (a < 0 || a >= ori_grid.i || b < 0 || b >= ori_grid.j || c < 0 || c >= ori_grid.k)
                            continue;
                        root_id = ori_grid.get_hex_id(a, b, c);
                        vector<double> temp_values;
                        vector<double> temp_coord;
                        vector<double> a_temp_values;
                        //获取root_coord中坐标点
                        temp_coord.push_back(ori_grid.hexdata[root_id].center.x);
                        temp_coord.push_back(ori_grid.hexdata[root_id].center.y);
                        temp_coord.push_back(ori_grid.hexdata[root_id].center.z);
                        //cout << "root_coord:" << temp_coord[0] << "," << temp_coord[1] << "," << temp_coord[2] << endl;
                        root_coord.push_back(temp_coord);
                        //获取root_coord中坐标点对应的杨氏模量、泊松比、密度、孔隙度s
                        a_temp_values.push_back(ori_grid.hexdata[root_id].poissonr);
                        a_temp_values.push_back(ori_grid.hexdata[root_id].youngmode);
                        a_temp_values.push_back(ori_grid.hexdata[root_id].rockden);
                        a_temp_values.push_back(ori_grid.hexdata[root_id].poro);
                        a_values.push_back(a_temp_values);

                        root_values[0].push_back(ori_grid.hexdata[root_id].poissonr);
                        root_values[1].push_back(ori_grid.hexdata[root_id].youngmode);
                        root_values[2].push_back(ori_grid.hexdata[root_id].rockden);
                        root_values[3].push_back(ori_grid.hexdata[root_id].poro);

                    }
                }
            }
            //获取泊松比最小值
            double min_poissonr = getExtremum(root_values[0])[0];
            //获取泊松比最大值
            double max_poissonr = getExtremum(root_values[0])[1];
            root_values[0] = normalize(root_values[0]);
            //获取杨氏模量最小值
            double min_yangmodel = getExtremum(root_values[1])[0];
            //获取杨氏模量最大值
            double max_yangmodel = getExtremum(root_values[1])[1];
            //对杨氏模量进行归一化
            root_values[1] = normalize(root_values[1]);
            //获取岩石密度最大值
            double min_rockden = getExtremum(root_values[2])[0];
            double max_rockden = getExtremum(root_values[2])[1];
            root_values[2] = normalize(root_values[2]);
            //获取孔隙度最大值
            double min_poro = getExtremum(root_values[3])[0];
            double max_poro = getExtremum(root_values[3])[1];
            if (min_poro != max_poro)
            {
                root_values[3] = normalize(root_values[3]);
            }
            //获取子网格的id值
            vector<int> children_id = ext_grid.children_id_of(root_i, root_j, root_k);
            //声明子网格坐标集
            vector<vector<double>> children_coord;

            //声明子网格体积集
            vector<double> children_volume(8);

            //声明子网格的8个顶点坐标集
            vector<vector<double>> children_mpoint(8);
            //vector<double> children_temp_mpoint(3);
            vector<double> children_temp_coord(3);
            vector<double> rbf_fourvalues(4);

            children_volume.resize(children_id.size());
            children_mpoint.resize(children_id.size());
            //获取子网格id值对应的坐标
            for (size_t children_i = 0; children_i < children_id.size(); ++children_i) {

                children_temp_coord[0] = ext_grid.hexdata[children_id[children_i]].center.x;
                children_temp_coord[1] = ext_grid.hexdata[children_id[children_i]].center.y;
                children_temp_coord[2] = ext_grid.hexdata[children_id[children_i]].center.z;
                //rbf_fourvalues = rbfInterpolateFourValues(root_coord, root_values, children_temp_coord);
                //ext_data[children_id[children_i]].poissonr = rbf_fourvalues[0];
                //ext_data[children_id[children_i]].youngmode = rbf_fourvalues[1];
                //ext_data[children_id[children_i]].rockden = rbf_fourvalues[2];
                //ext_data[children_id[children_i]].poro = rbf_fourvalues[3];

                children_coord.push_back(children_temp_coord);

                for (int childrenMpoint_i = 0; childrenMpoint_i < 8; ++childrenMpoint_i) {
                    children_mpoint[childrenMpoint_i].resize(3);
                    children_mpoint[childrenMpoint_i][0] = ext_grid.pdata[ext_grid.hexdata[children_id[children_i]].mpoint[childrenMpoint_i]].x;
                    children_mpoint[childrenMpoint_i][1]= ext_grid.pdata[ext_grid.hexdata[children_id[children_i]].mpoint[childrenMpoint_i]].y;
                    children_mpoint[childrenMpoint_i][2]= ext_grid.pdata[ext_grid.hexdata[children_id[children_i]].mpoint[childrenMpoint_i]].z;
                    //cout << "x: " << ext_grid.pdata[ext_grid.hexdata[children_id[children_i]].mpoint[childrenMpoint_i]].x << ", y: " << ext_grid.pdata[ext_grid.hexdata[children_id[children_i]].mpoint[childrenMpoint_i]].y << ", z: " << ext_grid.pdata[ext_grid.hexdata[children_id[children_i]].mpoint[childrenMpoint_i]].z << endl;
                    //children_mpoint[childrenMpoint_i].push_back(children_temp_mpoint);
                }
                //计算子网格体积
                //double children_temp_volume = hexahedronVolume(children_mpoint);
                children_volume[children_i] = hexahedronVolume(children_mpoint);
            }
            //声明原始网格8个顶点坐标集
            vector<vector<double>> root_mpoint;
            vector<double> root_temp_mpoint(3);
            //计算当前原始网格体积

            for (int rootVolume_i = 0; rootVolume_i < 8; rootVolume_i++) {
                root_temp_mpoint[0] = ori_grid.pdata[ori_grid.hexdata[i].mpoint[rootVolume_i]].x;
                root_temp_mpoint[1] = ori_grid.pdata[ori_grid.hexdata[i].mpoint[rootVolume_i]].y;
                root_temp_mpoint[2]= ori_grid.pdata[ori_grid.hexdata[i].mpoint[rootVolume_i]].z;
                root_mpoint.push_back(root_temp_mpoint);
            }
            double root_volume = hexahedronVolume(root_mpoint);

            //插值

            auto reservior_values = rbfExtrapolate(root_coord, root_values, children_coord, thread_num);
            reservior_values[0] = denormalize(reservior_values[0], min_poissonr, max_poissonr);
            //杨氏模量逆归一化
            reservior_values[1] = denormalize(reservior_values[1], min_yangmodel, max_yangmodel);
            //密度逆归一化
            reservior_values[2] = denormalize(reservior_values[2], min_rockden, max_rockden);
            if (min_poro != max_poro) {
                reservior_values[3] = denormalize(reservior_values[3], min_poro, max_poro);
            }
            //reservior_values[1] = destandardize(reservior_values[1], mean, stdev);
            //获取子网格中的密度和孔隙度计算守恒系数
            vector<double> children_rockden;
            double children_temp_rockden;
            vector<double> children_poro;
            double children_temp_poro;
           // cout << "reservior.size: " << reservior_values[0][0] << ", " << reservior_values[0][1] << ", " << reservior_values[0][2] << ", " << reservior_values[0][3] << ", children_id.size: " << children_id.size() << endl;


            for (size_t children_i = 0; children_i < children_id.size(); ++children_i) {

                //将插值结果赋值到子网格中
                ext_data[children_id[children_i]].poissonr = reservior_values[0][children_i];
                ext_data[children_id[children_i]].youngmode = reservior_values[1][children_i];

                if (ext_data[children_id[children_i]].n_child) {
                    for (int ii = 0;ii < ext_data[children_id[children_i]].n_child;ii++) {
                        ext_data[children_id[children_i]].child[ii].poissonr = reservior_values[0][children_i];;
                        ext_data[children_id[children_i]].child[ii].youngmode = reservior_values[1][children_i];;
                    }
                }
                
                ext_data[children_id[children_i]].rockden = reservior_values[2][children_i];
                ext_data[children_id[children_i]].poro = reservior_values[3][children_i];
                //保存密度和孔隙度的值
                children_temp_rockden = ext_data[children_id[children_i]].rockden;
                children_rockden.push_back(children_temp_rockden);
                children_temp_poro = ext_data[children_id[children_i]].poro;
                children_poro.push_back(children_temp_poro);
                vector<double> a_coord(3);
                if (ext_data[children_id[children_i]].poro > max_rbf_poro) {
                   max_rbf_poro = ext_data[children_id[children_i]].poro;
                }

            }

            //计算守恒系数
            vector<double> rockden_poro_coefficient = calculatedConservationCoefficient(root_volume, ori_grid.hexdata[i].rockden, ori_grid.hexdata[i].poro, children_volume, children_rockden, children_poro);
            //实现守恒
            for (size_t children_i = 0; children_i < children_id.size(); ++children_i) {
                //cout << "interpolation_rocken: " << ext_data[children_id[children_i]].rockden << endl;
                //cout << "rockden_poro_coefficient[1]: " << rockden_poro_coefficient[1] << endl;

                ext_data[children_id[children_i]].rockden = ext_data[children_id[children_i]].rockden * rockden_poro_coefficient[0];

                ext_data[children_id[children_i]].poro = ext_data[children_id[children_i]].poro * rockden_poro_coefficient[1];
                
                if (ext_data[children_id[children_i]].n_child) {
                    for (int ii = 0;ii < ext_data[children_id[children_i]].n_child;ii++) {
                        ext_data[children_id[children_i]].child[ii].rockden = ext_data[children_id[children_i]].rockden;
                        ext_data[children_id[children_i]].child[ii].poro = ext_data[children_id[children_i]].poro;
                    }
                }
                //cout << "conservation_rocken: " << ext_data[children_id[children_i]].rockden << endl;
                //cout << "conservation_poro: " << ext_data[children_id[children_i]].poro << endl;
            }
        }
    }

    void ExtrapolateInterpolate(CPgrid& ori_grid, finiteGrid& ext_grid, int thread_num)
    {
        //获取原始网格（Reservior）区域网格单元信息
        auto data = ori_grid.hexdata;
        //获取延拓后网格单元信息
        auto& ext_data = ext_grid.hexdata;
        //计算剖分后的网格单元中心坐标
#pragma omp parallel for schedule(static) num_threads(thread_num)
        for (int i = 0; i < ext_grid.hexsize; ++i)
        {
            ext_data[i].cal_center(&ext_grid.pdata[0]);
            // ext_data[i].youngmode = 200000;
        }
        //计算延拓后网格Reservior区域网格单元数量
        auto cell_reservior = ext_grid.topo_i * ext_grid.topo_j * ext_grid.topo_k;
        //插值矩阵构建选点
        PickingSpots pickingSpots(0.6);
        pickingSpots.setMin_Rate(0.001);

        //对储层上下两侧的延拓区域进行属性外推
#pragma omp parallel for schedule(static) num_threads(thread_num)
        for (int i = ext_grid.reservior_relative_i_start; i < ext_grid.reservior_relative_i_end; ++i)
        {
            auto selection = pickingSpots.Picking(ext_grid, i, "column");

            //获取selection坐标集
            vector<vector<double>> roofCoord(selection.size(),vector<double>(3,0.0));
            for (int m = 0; m < selection.size(); m++)
            {
                roofCoord[m][0] = ext_data[selection[m]].center.x;
                roofCoord[m][1] = ext_data[selection[m]].center.y;
                roofCoord[m][2] = ext_data[selection[m]].center.z;
            }
            //获取selection属性集
            vector<vector<double>> values(4,vector<double>(selection.size(),0.0));
            for (int i = 0; i < selection.size(); i++)
            {
                values[0][i] = ext_data[selection[i]].poissonr;
                values[1][i] = ext_data[selection[i]].youngmode;
                values[2][i] = ext_data[selection[i]].rockden;
                values[3][i] = ext_data[selection[i]].poro;
            }
            //属性归一化
            vector<vector<double>> normalized_values(4);
            for (int i = 0; i < values.size(); ++i)
            {
                normalized_values[i] = normalize(values[i]);
            }
            //获取目标区域单元数量
            // auto cells = (ext_grid.k_withext - ext_grid.topo_k) * ext_grid.j_withext;
            auto cells = ext_grid.k_withext * ext_grid.j_withext - ext_grid.topo_j * ext_grid.topo_k;
            //遍历延拓区域，获取延拓插值区域信息
            vector<vector<double>> targetPoints(cells, vector<double>(3));
            vector<int> target_index(cells);
            int cell_index = 0;
            int sum_ = 0;
            for (int ext_k = 0; ext_k < ext_grid.k_withext; ++ext_k)
            {
                if (ext_k < ext_grid.reservior_relative_k_start || ext_k >= ext_grid.reservior_relative_k_end)
                {
                    for (int ext_j = ext_grid.reservior_relative_j_start; ext_j < ext_grid.reservior_relative_j_end; ext_j++)
                    {
                        cell_index = ext_grid.get_hex_id(i,ext_j,ext_k, ext_grid.i_withext, ext_grid.j_withext, ext_grid.k_withext);
                        targetPoints[sum_][0] = ext_data[cell_index].center.x;
                        targetPoints[sum_][1] = ext_data[cell_index].center.y;
                        targetPoints[sum_][2] = ext_data[cell_index].center.z;
                        target_index[sum_] = cell_index;
                        ++sum_;
                    }

                }
                else
                    continue;

            }
            //属性外推计算
            // auto resualt_values = rbfExtrapolate(roofCoord, normalized_values, targetPoints, thread_num);
            auto resualt_values = rbfExtrapolate(roofCoord, normalized_values, targetPoints, thread_num);
            //属性逆归一化，还原属性到标准尺度
            for (int i = 0; i < resualt_values.size(); ++i)
            {
                denormalize(resualt_values[i],values[i]);
            }

            //延拓插值数据写回到网格中
            for (int cell = 0; cell < cells; ++cell)
            {
                ext_data[target_index[cell]].poissonr = resualt_values[0][cell];
                ext_data[target_index[cell]].youngmode = resualt_values[1][cell];
                ext_data[target_index[cell]].rockden = resualt_values[2][cell];
                ext_data[target_index[cell]].poro = resualt_values[3][cell];
            }
        }

        //计算剩余延拓区域属性
        pickingSpots.setCoarseFactor(0.4);
#pragma omp parallel for schedule(static) num_threads(thread_num)
        for (int k = 0; k < ext_grid.k_withext; ++k)
        {
            auto selection = pickingSpots.Picking(ext_grid, k, "row");
            //获取selection坐标集
            vector<vector<double>> roofCoord(selection.size(),vector<double>(3,0.0));
            for (int i = 0; i < selection.size(); i++)
            {
                roofCoord[i][0] = ext_data[selection[i]].center.x;
                roofCoord[i][1] = ext_data[selection[i]].center.y;
                roofCoord[i][2] = ext_data[selection[i]].center.z;
            }
            //获取selection属性集
            vector<vector<double>> values(4,vector<double>(selection.size(),0.0));
            for (int i = 0; i < selection.size(); i++)
            {
                values[0][i] = ext_data[selection[i]].poissonr;
                values[1][i] = ext_data[selection[i]].youngmode;
                values[2][i] = ext_data[selection[i]].rockden;
                values[3][i] = ext_data[selection[i]].poro;
            }
            //属性归一化
            vector<vector<double>> normalized_values(4);
            for (int i = 0; i < values.size(); ++i)
            {
                normalized_values[i] = normalize(values[i]);
            }
            //获取目标网格单元数量
            auto cells = ext_grid.i_withext * ext_grid.j_withext - ext_grid.topo_i * ext_grid.topo_j;
            //遍历延拓区域，获取延拓插值区域信息
            vector<vector<double>> targetPoints(cells, vector<double>(3));
            vector<int> target_index(cells);
            int cell_index = 0;
            int sum_ = 0;

            for (int i = 0; i < ext_grid.i_withext; ++i)
            {
                for (int j = 0; j < ext_grid.j_withext; ++j)
                {
                    //排除储层区域
                    if (i >= ext_grid.reservior_relative_i_start && i < ext_grid.reservior_relative_i_end && j == ext_grid.reservior_relative_j_start)
                        j = ext_grid.reservior_relative_j_end;
                    cell_index = ext_grid.get_hex_id(i,j,k, ext_grid.i_withext, ext_grid.j_withext, ext_grid.k_withext);
                    targetPoints[sum_][0] = ext_data[cell_index].center.x;
                    targetPoints[sum_][1] = ext_data[cell_index].center.y;
                    targetPoints[sum_][2] = ext_data[cell_index].center.z;
                    target_index[sum_] = cell_index;
                    ++sum_;

                }
            }


            //属性外推计算
            auto resualt_values = rbfExtrapolate(roofCoord, normalized_values, targetPoints, thread_num);
            //属性逆归一化，还原属性到标准尺度
            for (int i = 0; i < resualt_values.size(); ++i)
            {
                denormalize(resualt_values[i],values[i]);
            }

            //延拓插值数据写回到网格中
            for (int cell = 0; cell < cells; ++cell)
            {
                ext_data[target_index[cell]].poissonr = resualt_values[0][cell];
                ext_data[target_index[cell]].youngmode = resualt_values[1][cell];
                ext_data[target_index[cell]].rockden = resualt_values[2][cell];
                ext_data[target_index[cell]].poro = resualt_values[3][cell];
            }
        }

        cout << "rbfInterpolate end. " << endl;

    }

    // 克里金插值
    vector<vector<double>> krigingInterpolation(const vector<vector<double>>& roofCoord,
                                                const vector<vector<double>>& values,
                                                const vector<vector<double>>& targetPoints_ext,
                                                int thread)
    {
        cout << "Extrapolation of extended regional attributes start." << endl;
        cout << "The size of matrix A is " << roofCoord.size() << endl;
        cout << "The size of targetPoints_ext is " << targetPoints_ext.size() << endl;
        //插值矩阵大小：matrix_size * matrix_size
        size_t matrix_size = roofCoord.size();
        size_t value_types = values.size();
        size_t target_size = targetPoints_ext.size();

        double nugget = 0.0; // 块金效应
        double sill = 1.0;   // 基台值
        double range = 1.0;  // 范围


        //插值矩阵A
        DenseMatrix A(matrix_size,matrix_size,0.0);
        //值向量B
        vector<DenseMatrix> B(1, DenseMatrix(matrix_size, 1, 0.0));
        vector<vector<double>> result(value_types,vector<double>(target_size,0.0));

        // 构建插值矩阵A
        clock_t time_stt = clock();
#pragma omp parallel for schedule(static) num_threads(thread)
        for (int i = 0; i < matrix_size; ++i) {
            for (int j = 0; j < matrix_size; ++j) {
                if (i != j) {
                    double r = euclideanDistance(roofCoord[i], roofCoord[j]);
                    A.setValue(i, j, -semivariance(r, nugget, sill, range));
                }
            }
            double sum = 0.0;
            for (int j = 0; j < matrix_size; ++j) {
                if (i != j) {
                    sum += semivariance(euclideanDistance(roofCoord[i], roofCoord[j]), nugget, sill, range);
                }
            }
            A.setValue(i, i, sill + nugget + sum);
        }

        cout << "time use in Construct the interpolation matrix A is: " << (clock() - time_stt) / (double)CLOCKS_PER_SEC << "s" << endl;

        // 求解线性方程组 A * X = B
        DenseSolver solver(A,thread);
        time_stt = clock();
        for (size_t target_idx = 0; target_idx < target_size; ++target_idx) {
            DenseMatrix B(matrix_size, 1, 0.0);
            for (int i = 0; i < matrix_size; ++i) {
                double r = euclideanDistance(targetPoints_ext[target_idx], roofCoord[i]);
                B.setValue(i, 0, -semivariance(r, nugget, sill, range));
            }

            auto coefficients = solver.solve(B);
            cout << "time use in Solve AX = B is: " << (clock() - time_stt) / (double)CLOCKS_PER_SEC << "s" << endl;

            for (size_t value_idx = 0; value_idx < value_types; ++value_idx) {
                for (int i = 0; i < matrix_size; ++i) {
                    result[value_idx][target_idx] += coefficients.GetData_c(0)[i] * values[value_idx][i];
                }
            }
        }

        cout << "time use in get the interpolation results is       " << (clock() - time_stt) / (double)CLOCKS_PER_SEC << "s" << endl;
        cout << "rbfExtrapolate end." << endl;
        return result;
    }

    void assign(finiteGrid &grid,vector<vector<double>>values)
    {
        for (int k = 0; k < grid.k_withext; ++k) {
            for (int i = 0; i < grid.i_withext; i++) {
                for (int j = 0; j < grid.j_withext; j++) {
                    if (k < grid.reservior_relative_k_start) {
                        if (i >= grid.reservior_relative_i_start && i < grid.reservior_relative_i_end && j >= grid.reservior_relative_j_start && j < grid.reservior_relative_j_end) {
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poissonr = values[0][0];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].youngmode = values[0][1];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poro = values[0][2];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].rockden = values[0][3];
                        }
                        else {
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poissonr = values[1][0];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].youngmode = values[1][1];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poro = values[1][2];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].rockden = values[1][3];
                        }

                    }
                    else if (k < grid.reservior_relative_k_end) {
                        if (i >= grid.reservior_relative_i_start && i < grid.reservior_relative_i_end && j >= grid.reservior_relative_j_start && j < grid.reservior_relative_j_end) {
                            continue;
                        }
                        else {
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poissonr = values[2][0];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].youngmode = values[2][1];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poro = values[2][2];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].rockden = values[2][3];
                        }
                    }
                    else {
                        if (i >= grid.reservior_relative_i_start && i < grid.reservior_relative_i_end && j >= grid.reservior_relative_j_start && j < grid.reservior_relative_j_end) {
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poissonr = values[3][0];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].youngmode = values[3][1];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poro = values[3][2];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].rockden = values[3][3];
                        }
                        else {
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poissonr = values[4][0];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].youngmode = values[4][1];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].poro = values[4][2];
                            grid.hexdata[grid.get_hex_id(i, j, k, grid.i_withext, grid.j_withext, grid.k_withext)].rockden = values[4][3];
                        }
                    }
                }
            }
        }
    }


    //插值入口
    void Interp(CPgrid& ori_grid, finiteGrid& ext_grid, int thread_num)
    {
        rbfInterpolate(ori_grid, ext_grid, thread_num);
        cout << "rbfInterpolate end!" << endl;
        ExtrapolateInterpolate(ori_grid, ext_grid, thread_num);
        cout << "ExtrapolateInterpolate end!" << endl;
    }



}